;+
; NAME: 
;       IRDC_GAUSS
;
; PURPOSE: 
;       Use bivariate Gaussians (one each for BGPS source and IRDC)
;       built according to catalog fitted-ellipses to determine IRDC
;       assiciation probability.
;
; CATEGORY:
;       distance-omnibus IRDC association routine
;
; CALLING SEQUENCE:
;       prob = IRDC_GAUSS(struct [,\MAKE_PS])
;
; INPUTS:
;       STRUCT   -- BGPS structure containing coordinates, Bolocat
;                   fitted ellipse information, and others.
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       MAKE_PS  -- Create postscript files containing 3D rendering of
;                   the overlap of the two Gaussians.  Creates files
;                   in the './local/output' directory.
;
; OUTPUTS:
;       PROB     -- Probability of IRDC association based on the
;                   Overlap Integral between the two bivariate Gaussians.
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       IRDC_BLOCK -- Routine driver_irdc.pro reads in the Peretto &
;                     Fuller (2009) catalog & stores it in memory, as
;                     well as defines the bgps_irdc structure.
;
; MODIFICATION HISTORY:
;
;       Tue Oct 20 11:16:39 2009, Erik Rosolowsky <erosolo@A302357>
;
;		Written
;
;       Modified: 08/18/10, TPEB -- Changed from IRDC_LOOKUP to
;                                   IRDC_GAUSS, using bivariate
;                                   Gaussians for IRDC association.
;       Modified: 01/24/11, TPEB -- Added BGPS_IRDC structure element
;                                   for the peak tau tabulated in the
;                                   Peretto & Fuller (2009) catalog.
;       Modified: 02/01/11, TPEB -- Modified the routine to use the
;                                   irdc_full_data.sav IDL save file
;                                   (includes almost every column from
;                                   Peretto & Fuller (2009) catalog)
;                                   and also to create a Gaussian for
;                                   every IRDC within the grid size.
;       Modified: 02/15/11, TPEB -- Made routine compliant with new
;                                   version of IRDC_FULL_DATA.SAV that
;                                   tabulates IRDC sizes by semi-major axes.
;-

FUNCTION IRDC_GAUSS, s, MAKE_PS=make_ps
  
  ;; Call the common block containing the Peretto Catalog
  COMMON IRDC_BLOCK
  
  ;; Parse BGPS information into usable variables
  l = s.glon_cen       ;; Use centroid since making bivariate Gaussians
  b = s.glat_cen
  bmaj = s.maj
  bmin = s.min
  pa_bgps = s.pos_ang
  
  ;;===================================================================
  ;; NOTE: In irdc_full_data.sav, all IRDC coordinates and position  ;;
  ;; angles are in Galactic coordinates!                             ;;
  ;;===================================================================
  
  ;; Determine distances from BGPS source to all IRDCs in the catalog
  GCIRC, 2, l, b, pf_irdc.l, pf_irdc.b, dist
  min_dist = min(dist, hit)
  
  ;;===================================================================
  ;; As of 02/15/11, all ellipse sizes in DISTANCE-OMNIBUS are the   ;;
  ;;   semi-major and semi-minor axes!  This is to simplify things   ;;
  ;;   so I don't get confused.                                      ;;
  ;;===================================================================
  
  ;;===================================================================
  ;; Check for whether this is even worth it...                      ;;
  ;; ( dist < 3.0x sum of semi-majors )                                   ;;
  ;;   The value of 3.0x was chosen by investigating the amplitude   ;;
  ;;   of the overlap integral as a function of this value           ;;
  ;;   (P <~ 0.0001)                                                 ;;
  ;;===================================================================
  check_dist = 3.0 * (bmaj + pf_irdc[hit].maj)
  IF (min_dist GT check_dist) THEN RETURN,0b
  
  ;; Since this value is used to determine the grid size used for the
  ;;  calculation of the overlap integral, set it to be a minimum of 
  ;;  250" (the radius of the HVS postage stamp images).
  check_dist = check_dist > 250.
  
  
  ;;===================================================================
  ;; Make a grid every 2", centered on the IRDC coords (of some size);;
  ;; Have 2 copies of this grid (one for IRDC, one for BGPS).        ;;
  ;;===================================================================
  
  grid_space = 2.                                        ;; 2" grid spacing
  grid_size = long(round(2. * check_dist / grid_space))  ;; N grid points
  
  
  ;; Coordinate arrays for the grids
  l_grid  = (dindgen(grid_size)-round(grid_size/2.))*(grid_space/3600.) + $
            pf_irdc[hit].l
  b_grid  = (dindgen(grid_size)-round(grid_size/2.))*(grid_space/3600.) + $
            pf_irdc[hit].b
  
  
  
  ;;===================================================================
  ;; In order to select the most-associated IRDC, we consider each   ;;
  ;;   of the IRDCs within a 250" radius of the BGPS centroid.  So,  ;;
  ;;   we create a separate grid for each of these IRDCs and         ;;
  ;;   calculate the association amplitude, selecting the largest    ;;
  ;;   as the "real" association probability.  We do all this since  ;;
  ;;   the closest IRDC centroid may not be the one with the largest ;;
  ;;   major axis (and hence overlap the BGPS source).               ;;
  ;;===================================================================

  ;; Set parameters for the BGPS source & build the grid
  gauss_bgps = [l, b, bmaj, bmin, pa_bgps]
  bgps_grid = build_2d_gaussian(l_grid, b_grid, gauss_bgps)
  
  ;; Count the # of IRDCs within 250"
  ihit = WHERE( dist LE 250., nradius )

  ;; Catch errors if min(dist) > 250"
  IF nradius EQ 0 THEN BEGIN
     ihit = hit
     nradius = 1
  ENDIF
  
  ;; Array to hold the amplitudes for the various IRDCs
  amps = fltarr(nradius)
  
  ;; Loop through the IRDCs
  FOR kk=0, nradius-1 DO BEGIN
     
     ;; Build the IRDC grid
     gauss_irdc = [pf_irdc[ihit[kk]].l,pf_irdc[ihit[kk]].b,$
                   pf_irdc[ihit[kk]].maj,pf_irdc[ihit[kk]].min,$
                   pf_irdc[ihit[kk]].pa_irdc]
     irdc_grid  = build_2d_gaussian(l_grid, b_grid, gauss_irdc)
     
     ;;============================================
     ;; Setting up for the Overlap Integral as the
     ;;   statistic of choice for determining the
     ;;   association probability
     ;;============================================
     
     product = irdc_grid*bgps_grid
     
     amps[kk] = TOTAL(product)*TOTAL(product) / $
                ( TOTAL(irdc_grid*irdc_grid) * $
                  TOTAL(bgps_grid*bgps_grid) )
  ENDFOR
  
  
  ;; Select the proper bgps_irdc element
  ind = WHERE(bgps_irdc.cnum EQ s.cnum, nind)
  
  ;; Write the various items used to an IDL structure for use with the
  ;; Human Verification Slides
  bgps_irdc[ind].b_maj       = bmaj
  bgps_irdc[ind].b_min       = bmin
  bgps_irdc[ind].pa_b        = pa_bgps
  
  ;; Loop through the IRDCs indicated by nradius
  bgps_irdc[ind].n_irdc      = nradius
  bgps_irdc[ind].amplitude   = max(amps, imax)
  bgps_irdc[ind].hit_irdc    = imax
  FOR ll = 0L, nradius-1 DO BEGIN
     pfi = pf_irdc[ihit[ll]]
     bgps_irdc[ind].i_lcen[ll]   = pfi.l
     bgps_irdc[ind].i_bcen[ll]   = pfi.b
     bgps_irdc[ind].i_maj[ll]    = pfi.maj/2.
     bgps_irdc[ind].i_min[ll]    = pfi.min/2.
     bgps_irdc[ind].pa_i[ll]     = pfi.pa_irdc
     bgps_irdc[ind].taup[ll]     = pfi.tau_peak
     bgps_irdc[ind].contrast[ll] = 1. - (pfi.I_mir - pfi.I_min) / pfi.I_mir
     bgps_irdc[ind].amp[ll]      = amps[ll]
     bgps_irdc[ind].dist[ll]     = dist[ihit[ll]]
  ENDFOR     
  bgps_irdc[ind].majax_ratio = bgps_irdc[ind].i_maj[imax] / bmaj
    
  ;; If keyword set for creating postscript plots, do that now...
  
  IF KEYWORD_SET(make_ps) THEN BEGIN
     
     ;; Make masks for where each Gaussian is larger than the other
     bgps_mask = bgps_grid GT irdc_grid
     irdc_mask = irdc_grid GE bgps_grid
     
     plot_sum = (bgps_grid * bgps_mask) + (irdc_grid * irdc_mask)
     
     title=string(s.cnum,format="('BGPS #',I4)")
     
     set_plot,'ps'
     !p.font = 0
     mydevice,filename='local/output/irdc'+string(s.cnum,format="(I04)")+'.eps'
     
     shade_surf,plot_sum,l_grid,b_grid,title=title,$
                /xstyle,/ystyle,zrange=[0,1],/save,pixels=2048,$
                xtit='l (deg)',ytit='b (deg)'
     plots,pf_irdc[hit].l,pf_irdc[hit].b,1,psym=2,symsize=2,/t3d
     
     device,/close_file
     !p.font=-1
     set_plot,'x'
     
     print,title
     
  ENDIF
  
  ;; Return the Overlap Integral as the probability of association
  RETURN,bgps_irdc[ind].amplitude
  
END
